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Abstract. Generating accurate three dimensional planetary models is 
becoming increasingly important as NASA plans manned missions to re- 
turn to the Moon in the next decade. This paper describes a 3D surface 
reconstruction system called the Ames Stereo Pipeline that is designed to 
produce such models automatically by processing orbital stereo imagery. 
We discuss two important core aspects of this system: (1) refinement of 
satellite station positions and pose estimates through least squares bun- 
dle adjustment; and (2) a stochastic plane fitting algorithm that general- 
izes the Lucas-Kanade method for optimal matching between stereo pair 
images.. These techniques allow us to automatically produce seamless, 
highly accurate digital elevation models from multiple stereo image pairs 
while significantly reducing the influence of image noise. Our technique 
is demonstrated on a set of 71 high resolution scanned images from the 
Apollo 15 mission. 


1 Introduction 

Accurate, high resolution Lunar 3D maps will play a central role in NASA’s 
future manned and unmanned missions to the moon. These maps support land- 
ing site selection and analysis, lunar landing simulation & training efforts, and 
computer assisted landing systems. Furthermore, 3D digital elevation models 
(DEMs) provide valuable information to scientists and geologists studying lunar 
morphology. 

Several recent recent lunar satellite missions, including NASA’s Lunar Re- 
connaissance Orbiter, have returned stereo pairs with unparalleled resolution 
and image quality. However, historical data collected during the Apollo era still 
provide some of the best lunar imagery available today [1]. In fact, the Apollo 
Metric Camera system collected roughly 8,000 images covering roughly 20% of 
the lunar equatorial zone at a resolution of 10-m/pixel (Figure 1). The exten- 
sive coverage and relatively high resolution of this camera makes this data set 
extremely relevant in modern lunar data processing. 

In this paper, we introduce the Ames Stereo Pipeline, a C++ software frame- 
work for automated stereogrammetric processing of NASA imagery. We begin 
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Fig. 1 . Adjacent Apollo Metric Camera frames (e.g. AS15-M-1135 and AS15-M-1136 
shown here) overlap by 80%. This combined with the relatively wide field of view of 
the camera (74 degrees) results in ideal stereo angles between successive images. 


with an overview of this stereo reconstruction framework in Section 2. Then, spe- 
cific attention is given to two core components of the system: Section 3 describes 
the bundle adjustment approach for correcting extrinsic camera parameters and 
co-registering overlapping images; and Section 4 describes our sub-pixel accu- 
rate stereo correlation technique. Finally, in Section 5 we present the results of 
processing Apollo Metric Camera imagery. 

2 The Ames Stereo Pipeline 

The entire stereo correlation process, from raw input images to a point cloud or 
DEM, can be viewed as a multistage pipeline as depicted in Figure 2. 

The process begins with least squares Bundle Adjustment, which is described 
in Section 3, below. This produces corrected extrinsic camera parameters that 
are utilized by various camera modeling steps. 

Then, the left and right images are aligned using interest points or geometric 
constraints from the camera models. This step is often essential for performance 
because it ensures that the disparity search space is bounded to a known area. 
Next, a prepossessing filter such as the Sign of the Laplacian of the Gaussian filter 
is used, which has the effect of producing images that are somewhat invariant 
to differences in lighting conditions [2]. 

Following these pre-processing steps, we compute the disparity space image 
DSI(i,j,d x ,d y ) that stores the matching cost between a left image block cen- 
tered around pixel (i,j) and a right image block centered at position ( i — d x ,j — 
d y ). At this stage, the quality of the match is measured as the normalized cross 
correlation [3] between two 15x15 pixel image patches. We employ several opti- 
mizations to accelerate this computation: (1) a box filter-like accumulator that 
reduces duplicate operations in the calculation of DSI [4]; (2) a coarse-to-fine 
pyramid based approach where disparities are estimated using low resolution 
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Fig. 2. Flow of data through the Ames Stereo Pipeline 


images, and then successively refined at higher resolutions; and (3) partitioning 
of the disparity search space into rectangular sub-regions with similar values of 
disparity determined in the previous lower resolution level of the pyramid [4]. 

The DSI estimate just described efficiently computes integer estimates of 
disparity between the two images. These estimates are subsequently refined to 
sub-pixel accuracy using the technique described in Section 4. Finally, in con- 
junction with the bundle adjusted camera models, the sub-pixel disparity esti- 
mates are used to triangulate the location of 3D points as the closest point of 
intersection of two forward-projected rays emanating from the centers of the two 
cameras through the matched pixels. 

3 Bundle Adjustment 

The Apollo-era satellite tracking network was highly inaccurate by today’s stan- 
dards with errors estimated to be 2.04-km for satellite station positions and 
0.002 degrees for pose estimates in a typical Apollo 15 image [5]. Such errors 
propagate through the stereo triangulation process, resulting in systematic po- 
sition errors and distortions in the resulting DEMs (see Figure 3). These errors 
can be corrected using least-squares bundle adjustment. 

In bundle adjustment the position and orientation of the camera are deter- 
mined jointly with the 3D position of a set of image tie-points points chosen in 
the overlapping regions between consecutive images. Tie-points are automati- 
cally extracted using the SURF robust feature extraction algorithm [6] . Outliers 
are rejected using the RANSAC method and trimmed to 1000 matches that are 
spread evenly across the images. 

Our bundle adjustment approach follows the method described in [7] and 
determines the best camera parameters that minimize the projection error given 
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Fig. 3. Bundle adjustment is illustrated here using a color-mapped, hill-shaded DEM 
mosaic from Apollo 15 Orbit 33 imagery, (a) Prior to bundle adjustment, large disconti- 
nuities exist between overlapping DEMs. (b) After bundle adjustment, DEM alignment 
errors are no longer visible. 


by e = ^2j(Ik ~ I(Cj,Xi c)) 2 where I k are feature locations on the image 
plane, Cj are the camera parameters, and Xk are the 3D positions associated 
with features /&. I(Cj,Xk) is an image formation model (i.e. forward projection) 
for a given camera and 3D point. The optimization of the cost function uses the 
Levenberg-Marquartd algorithm. Speed is improved by using sparse methods 
described in [8]. 

To eliminate the gauge freedom inherent in this problem, we add two addition 
error metrics to this cost function to constrain the position and scale of the 
overall solution. First, e = JT — Cj ) 2 constrains camera parameters 
to stay relatively close to their initial values. Second, a small handful of 3D 
ground control points are chosen by hand and added to the error metric as 
e = ^2 k {Xj v cp — X k ) 2 to constrain these points to known locations in the lunar 
coordinate frame. In the cost functions discussed above, errors are weighted by 
the inverse covariance of the measurement that gave rise to the constraint. 

4 Sub-pixel Stereo Correlation 

Apollo images are affected by two types of noise inherent to the scanning process: 
(1) the presence of film grain and (2) dust & lint particles. The former gives 
rise to noise in the DEM values that wash out real features, and the latter 
causes incorrect matches or hard to detect blemishes in the DEM. Attenuating 
the effect of these scanning artifacts while simultaneously refining the integer 
disparity map to sub-pixel accuracy has become a critical goal of our system, 
and is necessary for processing real-world data sets such as the Apollo Metric 
Camera data. 

A common technique in sub-pixel refinement is to fit a parabola to the cor- 
relation cost surface in the 8-connected neighborhood around the integer dispar- 
ity estimate, and then use the parabola’s minimum as the sub-pixel disparity 
value. This method is easy to implement and fast to compute, but exhibits a 
problem known as pixel-locking: the sub-pixel disparities tend toward their inte- 


Lecture Notes in Computer Science 


5 


ger estimates and can create noticable ” stair steps” on surfaces that should be 
smooth [9], [10]. One way of attenuating the pixel- locking effect is through the 
use of a symmetric cost function [11] for matching the “left” and “right” image 
blocks. 

To avoid the high computational complexity of these methods another class 
of approaches based on the Lucas-Kanade algorithm [12] proposes an asymmet- 
ric score where the disparity map is computed using the best matching score 
between the left image block and an optimally affine transformed block from the 
right image. For example, the sub-pixel refinement developed by Stein et. al. [9] 
lets J#(ra, n ) and /^(i, j) be two corresponding pixels in the right and left image 
respectively, where i = m + d x , j = n + d y and d x , d y are the integer dispari- 
ties. They develop a linear approximation based on the Taylor Series expansion 
around pixel (i,j) in the left image 

I L (i + S X ,j + 5y) » I L (i, j ) + S x ^(i,j) + Sy^(i, j) (1) 

x 

where S x and S y are the local sub-pixel displacements. Let e(x,y) = In(x,y) — 
I L (i + S x , j + S y ) and W be an image window centered around pixel (m, n). The 
local displacements are not constant accross W and they vary according to: 

S x (hj) = a ii + hj + ci 

S y (iJ) = a 2 i + b 2 j + c 2 . (2) 

The goal is to find the parameters ai, &i, ci, a 2 , b 2l c 2 that minimize the cost 
function 

E (m,n)= ^2 .'/)»'(•'•• .'/)) 2 (3) 

(■ x,y)eW 

where w(x,y ) are a set of weights used to reject outliers. Note that the local 
displacements 5 x (i,j) and 5 y (i,j) depend on the pixel positions within the win- 
dow W. In fact, the values ai f bi,ci,a 2 ,b 2 ,c 2 that minimize E can be seen as 
the parameters of an affine transformation that best transforms the right image 
window to match the reference (left) image window. 

The shortcoming of this method is directly related to the cost function that it 
is minimizing, which has a low tolerance to noise. Noise present in the image will 
easily dominate the result of the squared error function, giving rise to erroneous 
disparity information. Recently, several statistical approaches (e.g. [13]) have 

emerged to show how stochastic models can be used to attenuate the effects 
of noise. Our sub-pixel refinement technique [14] adopts some of these ideas, 
generalizing the earlier work by Stein et. al. [9] to a Bayesian framework that 
models both the data and image noise. 

In our approach the probability of a pixel in the right image is given by the 
following Bayesian model: 

P(I R (m,n))= H M(I R (m,n)\I L (i + S x ,j + S y ),-^=)P(z = 0)+ (4) 

(x, y )ew ™ 

+ Af(I R (m,n)\fj, n ,(T n )P(z = 1) 
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The first mixture component (z = 0) is a normal density function with mean 
I L (i + j + S v ) and variance Zf— : 

y y/9xy 

P(I R (m,n)\z = 0 ) = Af(I R (m,n)\I L (i + S X , j +<5 y ), ~^=) (5) 

v 9xy 

The i factor in the variance of this component has the effect of a Gaussian 

v 9xy 

smoothing window over the patch. With this term in place, we are no longer 
looking for a single variance over the whole patch; instead we are assuming the 
variance increases with distance away from the center according to the inverted 
Gaussian, and are attempting to fit a global scale, a p . This provides formal 
justification for the standard Gaussian windowing kernel. 

The second mixture component (z = 1) in Equation 5 models the image noise 
using a normal density function with mean fjL n and variance cr n : 

P(I R (m,n)\z = 1) = Af(I R (m, rc)|/x n ,<j n ) (6) 

Let I R (m,ri) be a vector of all pixels values in a window W centered in pixel 
(m, n) in the right image. Then, 

P(I fl (m,n))= J] P(I R (x, y)) (7) 

(x,y)£W 

The parameters A = {ai, &i, ci, a 2 , 62 , C 2 , cr p , /i n , <j n } that maximize the model 
likelihood in Equation 7 are determined using the Expectation Maximization 
(EM) algorithm. Maximizing the model likelihood in Equation 7 is equivalent to 
maximizing the auxiliary function: 

Q(0) = E A ‘) k » ^l A ) 

k 

= EE P(k\I R (x, y), At) log P(I R (x, y)\k, X)P(k\X) ( 8 ) 

k x,y 

Note that the M step calculations are similar to the equation used to deter- 
mine the parameters ai, &i, ci, a 2 , 62 , in the method presented in [9], except 
here the fixed set of weights is replaced by the a posteriori probabilities computed 
in the E step. In this way, our approach can be seen as a generalization of the 
Lucas-Kanade method. The complete algorithm is summarized in the following 
steps: 

— Step 1 : Compute ^(i, j), ^(i, j) and the I R (x,y) values using bilinear 
interpolation. Initialize the model parameters A. 

— Step 2: Compute iteratively the model parameters A using the EM algorithm 
(see [14] for details). 

— Step 3: Compute S x (i,j) and 5 y (i,j) using Equation 2. 

— Step 4: Compute a new point = (x,y) + (S x ,5 y ) and the I R (x f ,y') 

values using bilinear interpolation. 
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Fig. 4. Hadley Rille and the Apollo 15 landing site derived from Apollo Metric Camera 
frames AS15-M-1135 and AS15-M-1136. (a) superimposed over the USGS Clementine 
base map, (b) oblique view. 


— Step 5: If the norm of (S x ,5 y ) vector falls below a fixed threshold the iter- 
ations converged. Otherwise, go to step 1. 

Like the computation of the integer DSI , we adopt a multi-scale approach 
for sub-pixel refinement. At each level of the pyramid, the algorithm is initialized 
with the disparity determined in the previous lower resolution level of the pyra- 
mid. This allows the subpixel algorithm to shift the results of the integer DSI 
by many pixel if a better match can be found using the affine, noise-adapted 
window. 

5 Results 

The 3D surface reconstruction system described in this paper was tested by 
processing 71 Apollo Metric Camera images from Apollo 15. Specifically, we 
chose frames from orbit 33 of the mission, which includes highly overlapping 
images that span approximately 90 degrees of longitude in the lunar equatorial 
region. This exercised our algorithms across a wide range of different terrain 
and lighting conditions. Figure 4 shows the final results in the vicinity of Hadley 
Rille: the Apollo 15 landing site. 

Tests were carried out on a 2.8-GHz, 8-core workstation with 8-GB of RAM. 
Stereo reconstruction for all 71 stereo pairs took 2.5 days. In the end, the results 
were merged into a DEM at 40-m/pixel that contained 73,000 x 20,000 pixels. 

5.1 Bundle Adjustment 

Bundle adjustment was carried out as described in Section 3. Initial errors and 
results after one round of adjustment are shown in columns two and three of 
Table 1, respectively. Subsequently, any tie-point measurements with image- 
plane residual errors that were greater than 2 standard deviations from the 
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mean residual error were thrown out. Bundle adjustment was run a second time 
yielding slightly improved results shown in column four of the table. 

To constrain the scale and absolute position of the solution, 7 ground control 
points were selected in a triangle wave pattern across the extent of the orbit to tie 
specific image pixels to known positions in the lunar coordinate frame. The sigma 
weights for ground control points were based on the resolution of the underlying 
base map from which ground control points were derived. These were 300-m on 
the surface and 500-m normal to the surface. Furthermore, the camera station 
position and pose estimates were constrained to stay close to their initial values 
based on radio tracking data. Sigma weights for camera parameters were 2-km 
for position, and 0.01 radians for pose. These values were drawn from historical 
estimates of Apollo tracking network accuracy as previously discussed. 


Residual Reconstruction 

Initial 

After Round 1 

After Round 2 

Image Plane 

0.444-mm 

0.012-mm 

0.0075-mm 

Camera Position 

0-km 

1.31-km 

1.31-km 

Camera Orientation 

0-mrad 

9.0-rad 

9.1-mrad 

Ground Control Point 

0-m 

481-m 

465-m 

Triangulation Error 

911-m 

24.1-m 

15.58-m 


Table 1 . Residual error at various stages of bundle adjustment. Residual error in the 
image plane decreases as image tie-point constraints are satisfied. This improvement is 
made possible as residual “error” for camera position, orientation and ground control 
points increase to compensate. Triangulation is a measure of the average distance 
between the closest point of intersection of two forward projected rays for a set of 
tie-points. Its decrease indicates a substantial improvement in the self-consistency of 
the DEMs in the data set. 


5.2 Subpixel Correlation 


Film grain and the dust particles are inherent to the scanning process and can 
significantly limit the accuracy of the stereo processing system. One example 
where dust particle noise occurs in one of the stereo pair images is shown in 
detail in Figure 5 (a) and (b). Figure 5 (c) illustrates the integer disparity map 
obtained by running the fast discrete correlation method described in Section 2. 
Figure 5 (d), and (e) compares the horizontal sub-pixel disparity maps obtained 
using the parabola method and the Lucas-Kanade method with the Bayesian 
approach we introduced in Section 4. The Bayesian approach reduces the “stair- 
stepping” artifacts apparent in the results from the parabola method. It also 
demonstrates a degree of immunity to the noise introduced by the speck of dust 
in (a). 
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(d) 


(e) 


Fig. 5. (a) Left Image (with a speck of dust), (b) Right Image, (c) Horizontal integer 
disparity map, (d) Horizontal disparity map using the parabola method, (e) Horizontal 
disparity map using the Bayesian approach. 


6 Conclusions and future work 


This paper has introduced a novel statistical formulation for optimally determin- 
ing stereo correspondence with subpixel accuracy while simultaneously mitigat- 
ing the effects of image noise. Furthermore, we have successfully demonstrated 
a significant improvement to the geometric consistency of the results after us- 
ing least squares bundle adjustment. These techniques were successfully used 
to process a large, real-world corpus of images and were found to produce use- 
ful results. However, this was a preliminary demonstration of these capabilities, 
and much work remains to quantify residual errors and characterize the degree 
of noise immunity in our new correlation algorithm. Further research will be 
directed towards building a comprehensive error model for the 3D surface re- 
construction process that can be used to systematically test our system under a 
variety of different conditions and inputs. Ultimately, we hope to use our tech- 
nique to process the full collection of over 8,000 Apollo Metric Camera stereo 
pairs. 
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